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In light of the evidence that charge order coexists with d-wave superconductivity in the under- 
doped cuprate superconductors, we investigate the manner in which such charge order wiU influence 
the quasiparticle excitations of the system and, in particular, the low-temperature transport of heat 
by those quasiparticles. We consider a d-wave superconductor in which the superconductivity co- 
exists with charge density wave order of wave vector {■n/a,Q). While the nodes of the quasiparticle 
energy spectrum survive the onset of charge order, there exists a critical value of the charge density 
wave order parameter beyond which the quasiparticle spectrum becomes fully gapped. We perform 
a linear response Kubo formula calculation of thermal conductivity in the low temperature (uni- 
versal) limit. Results reveal the dependence of thermal transport on increasing charge order up to 
the critical value at which the quasiparticle spectrum becomes fully gapped and thermal conduc- 
tivity vanishes. In addition to numerical results, closed-form expressions are obtained in the clean 
limit for the special case of isotropic Dirac nodes. Signatures of the influence of charge order on 
low-temperature thermal transport are identified. 

PACS numbers: 74.25.Fy, 74.72.-h 



I. INTRODUCTION 

The low energy quasiparticle excitations of the d-wave 
superconducting phase of the high- Tc cuprate supercon- 
ductors are massless anisotropic Dirac fermions,^ These 
Dirac quasiparticles are easily excited in the vicinity of 
the four nodes, the four points on the two-dimensional 
Fermi surface where the superconducting order param- 
eter vanishes. The dominant carriers of heat at low 
temperature, quasiparticles are efficiently probed via low 
temperature thermal conductivity measurements, which 
have been performed extensively over the past decade. 
Theorj*2i^i^i^i^ii!^ has shown that the massless Dirac en- 
ergy spectrum yields a low temperature limit where ther- 
mal conductivity is remarkably independent of disorder 
for small impurity density. In this limit, known as the 
universal limit, thermal conductivity per Cu02 plane 
depends only on fundamental constants and the ratio 
of the Fermi velocity, vp: to the fc-space slope of the 
superconducting order parameter at the nodal points, 
VA- Experiments^iiaiiiJ^iiliMiiSiie^rzj^ l^^^^e demon- 
strated this disorder-independence and used this result 
to extract the anisotropy ratio, a = vf/va, from low- 
temperature thermal transport data. 

Over the past few years, there has been a signifi- 
cant effort to grow and measure high quality cuprate 
samples in the underdoped regime of the superconduct- 
ing phase, as well as the pseudogap phase that re- 
sults from underdoping even further. Several experi- 
mental groups have used high-resolution scanning tun- 
neling microscopyi2^'2i^iS3^'25^^^^^ to exam- 
ine the electronic states of the underdoped cuprates at 
the atomic scale. These experiments, amongst others^, 
have provided evidence that charge order coexists with 
d-wave superconductivity (dSC) in these materials. Fur- 



thermore, it has been shown theoreticall y^^'^'^i that co- 
existing charge order can significantly affect the quasi- 
particle spectrum of the superconductor, leading the sys- 
tem to become fully gapped for charge order of sufficient 
magnitude. If the quasiparticles are fully gapped (no 
nodes), the dominant carriers of heat at low temperature 
are frozen out, which should have a dramatic effect on 
the universal-limit thermal conductivity. The purpose of 
our current analysis is to study the nature of this effect. 

We consider a particularly simple form of charge or- 
der, a conventional s-wave charge density wave (CDW) 
with a fc-independent order parameter and a wave vector 
Q = (7r/a, 0) that doubles the unit cell. While the charge 
order in the underdoped cuprates may be of a more com- 
plex type, this simple model provides a place for us to 
start studying, phenomenologically, the effect of charge 
order on thermal transport in a d-wave superconductor. 
Furthermore, since the experimentally observed" CDW 
has a wave vector close to (7r/2a,0), it will generically 
have a second harmonic near (tt/o, 0). This harmonic 
can couple efficiently to the nodal quasiparticles because 
its wave vector nearly spans the separation between the 
nodes. "^^ The calculations presented in this paper can 
then be viewed as applying to this second harmonic. 

While the charge order in the cuprates may turn on 
with underdoping, we simply add a CDW term to the 
dSC Hamiltonian and turn on the charge order by hand, 
by increasing the magnitude of the CDW order parame- 
ter. We then calculate the universal limit thermal con- 
ductivity of the combined system, evaluating the effect 
of coexisting charge order on thermal transport. Our 
goal is to identify signatures of the onset of charge or- 
der which may be observed with underdoping in low- 
temperature thermal conductivity measurements of the 
underdoped cuprates. Evidence of the breakdown of 
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universal thermal conductivity at low doping, possibly 
due to the onset of charge order, has already been seen 
experimentally^i^'^i^, and our results might therefore 
shed light on these studies. 

We begin in Sec. [Til by writing down the com- 
bined Hamiltonian, calculating the resulting energy spec- 
trum, and discussing the charge-order-induced transi- 
tion whereby the spectrum can become fully gapped. In 
Sec. mil we calculate the Green's function and thermal 
current operator, define our model of disorder, and use 
a diagrammatic Kubo formula approach to obtain an in- 
tegral form for the thermal conductivity tensor. For the 
special case of a clean system (no disorder) with isotropic 
nodes {vp ~ wa) the remaining fc-space integration can 
be performed analytically. This case is considered in 
Sec. II VI where a closed-form solution is obtained for the 
thermal conductivity tensor as a function of the mag- 
nitude of the charge order. The more general case of 
nonzero disorder and anisotropic nodes is considered in 
Sec. |V] via a numerical computation, and the effect of 
disorder and nodal anisotropy is discussed. Conclusions 
are presented in Sec. IVB 



II. COEXISTING dSC AND CDW ORDER 

A. Hamiltonian 

Following Ref.m, we consider a model Hamiltonian for 
a d-wave superconductor with coexisting charge order: 



H — Hq + Hdsc + HcDw 
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FIG. 1: Charge order of wave vector Q = {n/a,0) doubles 
the unit cell and thereby halves the Brillouin zone. With in- 
creasing charge density wave order parameter, ■(/), the nodes of 
the energy spectrum, and their images in the second reduced 
Brillouin zone (shaded) , approach the reduced Brillouin zone 
edges (dotted), coUiding for ij^ = ij^c, beyond which the spec- 
trum is fully gapped. 



The charge density wave has the effect of doubling 
the unit cell and therefore halving the effective Brillouin 
zone, as shown in Fig. [TJ By defining a four-component 
extended-Nambu vector 



"fe+QT' '^-k-Ql 



(5) 



consisting of particle and hole operators at k and k-|- Q, 
we can express the Hamiltonian in a compact 4x4 matrix 
notation: 



ka 



(4) 



(6) 



Momenta are summed over the Brillouin zone of a two- 
dimensional square lattice of lattice constant a. Hq + 
Hdsc is the mean-field BCS Hamiltonian for electron 
dispersion and superconducting order parameter A^, 
which is taken to have d-wave symmetry (for example, 
Ak = Ao{coskxa — coskyo) /2). Hcuw denotes a charge 
density wave of wave vector Q with CDW order parame- 
ter ip. While it is possible to consider density- wave states 
of nonzero angular momentum**^ by taking ip to be com- 
plex and fc-dependent, we shall focus here on the effect 
of a conventional s-wave CDW corresponding to a site- 
centered charge modulation in the z-direction of wave- 
length twice the lattice constant. That is, we take ip to 
be a real. A:- independent parameter and set Q = (vr/a, 0). 



where 



Hk 



ei Ai V 

Ai -ei -tjj 

£2 A2 

-tl> A2 -£2 



(7) 



and subscript 1 denotes k and subscript 2 denotes k+ Q. 
The prime indicates that the momentum sum is restricted 
to the reduced Brillouin zone. Note that the upper-left 
and lower-right 2x2 blocks of Hk are simply the Nambu 
space Hamiltonian at k and k -I- Q respectively. The 
CDW order parameter couples these two sectors. 
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B. Energy Spectrum and Nodal Collision 

The energy spectrum of the fcrmionic excitations of 
this system of coexisting dSC and CDW order is obtained 
by solving for the (positive) eigenvalues of Hf^. Doing so, 
we find that 



2V/) 



(8) 



± 



{el + Al~el 



A 



1/2-1 1/2 



((ei+e2)2 + (Ai-A2)2)] } 



For ip = 0, these solutions reduce to the energy spectra of 
the quasiparticle excitations of th e d-wave superconduc- 

By construc- 



tor, E'l and -E^+g, where 



Al. 



tion, the quasiparticle energies drop to zero at four points 
in the Brillouin zone, the intersection of the Fermi surface 
with the lines = ±ky. These are the nodal points, or 
nodes, of the d-wave superconductor. To model the situ- 
ation in the cuprates, the nodes are taken to be a distance 
kp from the origin, inside of the (±7r/2a, ±7r/2a) points 
by a small distance fcg where /eq = Tr/v^a — kp <^ kp- 

As the charge density wave is turned on, the nodal 
structure of the excitation spectrum initially survives, 
since the CDW wavevector, Q, is not commensurate with 
the internodal distance. '^^■^^ With increasing ip, the nodes 
move toward the reduced Brillouin zone edge along the 
trajectory sketched in Fig. [T] Also plotted in this figure 
is the trajectory of the image of each node translated by 
Q into the second reduced Brillouin zone. At a critical 
value of the CDW order parameter, ip — ipc, the nodes 
collide at the reduced Brillouin zone edge and the energy 
spectrum becomes fully gapped. For tp > tpc-, the mini- 
mum values of the excitation spectra are nonzero. Hence, 
the nodes have vanished. 

To determine the points in /c-space at which this nodal 
collision occurs, we need only solve for the points at 
which a zero of Ek coincides with a reduced Brillouin 
zone edge. We define the tp for which this occurs to be 
ipc- For example, node #1 (located in the upper-right 
quadrant) will collide somewhere along the reduced zone 
boundary at k^ = n/2a. Setting E{kx = iT/2a,ky) — 
and noting that both and are even functions of 
kx^ we find that the collision point must satisfy ei = 
£2 = V'c and Al = A2 = 0. Near node #1, the lat- 
ter condition yields k^ — ky, so the collision point is 
kc — (7r/2a, 7r/2a). Equivalent arguments for each of the 
four quadrants reveal that the four collision points are 
located at (±7r/2a, ±7r/2a). Defining local coordinates 
ki and fc2 about each of the collision points, as shown in 
Fig. [2I we can write 



ei=vp{kQ + ki) Al = WAfc2 

^2 = VF{ko + k2) A2 = WAfcl 



(9) 



where vp is the Fermi velocity and va. is the slope of 
the gap at the node. Note that in writing these linear 
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FIG. 2: Local coordinates, k\ and ^2, defined about each of 
the four nodal collision points, kc = (±7r/2a, ±7r/2a). The fci- 
axes, perpendicular to the Fermi surface, define the direction 
of increasing electron dispersion, e^. The fc2-axes, parallel to 
the Fermi surface, define the direction of increasing supercon- 
ducting order parameter, Afc. 



relations, we have assumed that /cq is small enough that 
the spectrum of the d-wave superconductor is still linear 
in the vicinity of the collision points. At the collision 
points (fci = fc2 = 0), ei = e2 = wf^Oj which requires 
that tpc — vpko- Switching to scaled coordinates, pi = 
y/vpvXki and p2 = ^JvFV/^k2 , yields 



ei = + Al = p2 

£2 = V^c + \/ap2 A2 = pi 



(10) 



where a = vp/va- This notation provides a convenient 
framework with which to proceed with the thermal trans- 
port calculation. 



III. TRANSPORT CALCULATION 

Given the Hamiltonian defined by Eqs. ([7]) and ([TT 
we can calculate the thermal conductivity, and its de- 
pendence on the charge density wave order parameter, 
via Kubo formula. 



A. Green's Function 

We begin by computing the Matsubara Green's func- 
tion. In the extended-Nambu basis of Eq. the bare 
Green's function is a 4 x 4 matrix obtained through in- 
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version of the Hamiltonian 
It takes the form 



1 



Ga Gb 

Gc Gd 



Ga = ((zu;)2-e2- A2)[iu; + eir3 + Airi] 
-^l)'^[iuj - e2T3 + A2T1] 



(11) 



(12) 



(13) 



Gb = ij[iuj{ei + 62) + {{iuj)^ + £162 - - ip'^)T3 

+ (eiA2 + e2Ai)Ti - iuj{Ai - A2)(ir2)] (14) 

Gc ^P[iuj{ei + 62) + {{iujf + eie2 - A1A2 - ip^)T3 
+ (eiA2 + e2Ai)Ti + iw(Ai - A2)(iT2)] (15) 



Gd = {{iuf ~ el - Al)[iuj + e2T3 + A2T1] 
— eiT3 + AiTi] 



(16) 



Gdcn = {ei + A\+ij^-{iujf){el + Al+^^-{iuf) 
-^'((ei+e2)' + (Ai-A2)2) (17) 

where Gdon is a scalar and Ga, Gb, Gc, and Gd are 2x2 
matrices expressed in terms of particle-hole-space Pauli 
matrices, r^. 

In the presence of disorder, we must include the impu- 
rity contribution to the self-energy via Dyson's equation 



G-^ = G° 



(18) 



The self-energy, S, is a 4 x 4 matrix in the extended- 
Nambu basis, but for simplicity, we consider here only 
the scalar term 



(19) 



and postpone discussion of the effects of off-diagonal 
self-energy terms to a separate publication*^. Then the 
dressed Matsubara Green's function is simply 

G{k,iuj) = [{iuj - J:{iuj))l - Hkr^ ^ G°{k,iuj~J:{iuj)) 

(20) 

For our calculation of the zero-temperature thermal con- 
ductivity, we will require only the imaginary part of the 
zero-frequency retarded Green's function, ImG^(fc,tt' 
0). Continuing itu to + i6 and taking the w — *■ limit, 
the zero-frequency retarded self-energy is just a negative 
imaginary constant, — iFq, and we find that 



ImG^(fc,w^O) 



Gdc 



Gjj G^ 
Gc G J 



G„ 



-ro(r2 + v^2 



(21) 



(22) 



G, = ijTo [(ei + £2) - (Ai - A2)(iT2)] (23) 
g1' = ^pVo [(ei + £2) + (Ai - A2)(iT2)] (24) 
Gd^-Toirl+4'^ + 4 + Al) (25) 

Gdcn = iTl+^4j^ + el + Al)iTl+i,^ + el + Al) 

-V^'((ei+e2)' + (Ai- A2)2) (26) 

where Tq is the zero-frequency impurity scattering rate 
(the impurity-induced broadening of the spectral func- 
tion). 

B. Current Operator 

Next we must calculate the quasiparticle current oper- 
ator for this system of coexisting d-wave superconductor 
and charge order. We note that quasiparticles carry a 
well-defined heat and spin. Thus, where a quasiparticle 
goes, so goes its heat and spin. Though the quantity we 
require is the thermal current, we will proceed by cal- 
culating the spin current operator (which is technically 
simpler) obtaining the thermal current operator by cor- 
respondence. 

The spin current operator, j'', is obtained via continu- 
ity with the spin density operator, p". 



1 



(27) 



Fourier transforming and taking the zero-wavevector 
limit yields a recipe for calculating J^^q, which is the 
operator we will need for the transport calculation. 



(28) 



Defining and re-expressing the spin density operator in 
various forms, we note that 

Pq = ^ ^'^'^l' a'^k' +qa 
k' (J 
I 

k' 



(29) 



where Sa = ±s, s = 1/2, d^a = Ck+Qa, ^'fc is the four- 
component extended-Nambu vector defined in Eq. ([5]), 
and the prime restricts the wave vector sum to the re- 
duced Brillouin zone. In the same notation, the Hamil- 
tonian takes the form 
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we can calculate the thermal conductivity via the Kubo 
formula^ 



d-kicLki) 



d-kidk]) 
(30) 



KiT) 



liui „„ 

where the retarded current-current correlation function is 
obtained from the Matsubara function via analytic con- 



(38) 



Using fermion anticommutation relations to evaluate the 
commutator in Eq. we find that 



J5 



k 



VF fc+QT3 + VA k+QTi 



tinuation. 

u^iQ) ^fi^iin^n + id) (39) 

^fe+QIn what follows, we neglect vertex corrections, calculat- 



(31) 

where vpk = dek/dk, v^k = dAk/dk, and v^k = 
dip / dk. For the case we consider, ip is fc-independent, 
so v^fe is precisely zero and the spin current operator is 
block diagonal in the extended-Nambu basis. 

In the vicinity of each of the four collision points (the 
regions we will always be considering), wpk points along 
the locally-defined Aii-direction and va*: points along the 
locally-defined fc2-direction, as shown in Fig. [2] There- 
fore, shifting by wave vector Q = (tt/o, 0) from k to k-l-Q 
flips the sign of the x-component of each velocity while 
preserving the y-component. That is, the components 
satisfy = Tj^u^fe and ^Afc+Q = ^i«Afc for 



So we can write 



J5 



-1 for 
+1 for 



^ *i [VA/F + Va/a] 



I = X 

i = y 



k+Q 



where 



VAf A 



^Afe 



Ml 







Ml 



<km 







(32) 

(33) 

(34) 
(35) 

(36) 



Finally, we note that since the same quasiparticles that 
carry the spin also carry the heat, the thermal current 
operator, j", will have the same structure as the spin 
current operator. In the zero-wavevector, zero-frequency 
limit that we will require. 



Jo = 



C. Thermal Conductivity 



OJ + O 



(37) 



Given the Green's function, thermal current operator, 
and coordinate system defined in the previous sections. 



ing the bare bubble current-current correlation function 
using the Matsubara formalism^. It has been shown 
previously^ that vertex corrections are negligible for the 
d-wave superconductor case (without charge order) and 
the contribution of vertex corrections to the present case 
will be considered in the a separate paper^i. 

Evaluating the bare bubble Feynman diagram shown 
in Fig. [3] yields 



1 



if? ,o 

— ) 
2 ' 



X Tr [G{k, iLo)wMG{k, ilo -f ir2)vA/] 



(40) 



where va/ = va^f + va/ a is a vector in coordinate space 
and a matrix in extended-Nambu space, the Green's 
functions are dressed with disorder, the w-sum is over 
fcrmionic Matsubara frequencies, the fc-sum is restricted 
to the first reduced Brillouin zone, the trace is over 
extended-Nambu space, and j3 — l/ksT. We expand 
the fc-sum from the reduced Brillouin zone to the full 
(original) Brillouin zone, which double-counts and there- 
fore requires division by 2. Since the summand is sharply 
peaked in the vicinity of the four nodal collision points, 
we then replace the fc-sum by four integrals over local 
scaled coordinates, pi andp2, defined (in Sec lIIBj l about 
each of these points. 



1 



4 



P 



(2tt)'^vfva 



(41) 



Making use of a spectral representation of the matrix 
Green's function 



G(p,iuj) 



Eq. (|40|) becomes 



duji ■ 



-^ImG^(p,c^i) 
iuj — uJi 



(42) 



1 



J (27r 



P 



dujiduj2S{iD,) 



Tr 



^G;(p,c.i)vj^^) 



where 



s{in) = - 



1 



1 



(43) 



(44) 
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FIG. 3: Feynman diagram depicting the bare bubble ther- 
mal current-current correlation function, n„(iS7). The ther- 
mal current operator sits on each vertex and each propagator 
denotes a Green's function dressed with disorder self-energy. 



and vj(^^ is the value of vm in the vicinity of collision 
point j. (Note that while the spectral representation de- 
fined in Eq. (|42p is valid for the case of real ip that we 
are considering, it would not be valid if tp, and therefore 
Hk, was complex. The subtleties of this are discussed in 
detail in the Appendix.) 

Computing the Matsubara sum in Eq. (|44)l via contour 
integration (see Refs. andU for discussion of technical 
points), continuing + iS to obtain the retarded 

function, and taking the imaginary part, we find that 

S]^{Q) — 7r(tJi-|-— )^(ni7(tJi + f7)-ni?(uJi))(5(aJi+ri-aj2) 

(45) 

where np{x) = l/{e^^ + 1) is the Fermi function the 
double-prime indicates the imaginary part. Then taking 
the 57 ^ limit in Eq. (|38p yields an expression for the 
thermal conductivity tensor 



«(r) 

T 

where 



-1 



dco 



W \ 2 Qfip 

duj 



(f) 



47r 



TrR{p,uj) 



(46) 



i?(P,^) 



^Gfl,(p,c.i)v^^^G^(p,c.2)vi^) (47) 



and taking the T — > limit yields 

fp 
47r 



^0 _ fc| 



Tri?(p,0) 



(48) 



Here we have used the fact that, for low T, 
{u! /T)'^{—dnp /dco) is sharply peaked at w = and 



/•°° , /tJ\2 / dnF\ TT^kl 



(49) 



Noting that at each collision point, wp and va point 
along the local ki and directions respectively (as de- 
fined in Fig. [2]) and performing the sum over collision 
points in Eq. (|47p. we find that 

Ruij>, 0) = 2vlG"jiMlG"jiMl + 2vIg"^MIG"j^MI 
+2Tj,vpVAiG'RM;G'^Ml + GrMIGrMI) (50) 

for i = {x, y\ while 



i?,„(p,0) = i?j,,(p,0) = 



(51) 



where rji, M|, and Ml are defined in Eqs. ([32)l and (|36p . 
Plugging in the Green's function from Eqs. (PT] - [^5)) and 
taking the trace over the 4x4 extended-Nambu space 
yields 



T 



^00 

T 



(Pp Ni + r]iN2 



D 



(52) 



Ni ^2A[{A + B + el + Aff + {A + B + el + A^)^] 

(53) 

N2 - AAB [(ei + 62)2 - (Ai - As)^] (54) 



D = [{A + B + el + Al){A + B + el + A 
~ B((ei+e2)' + (Ai-A2)2)]' 



where 



Kqo _ k% /vp ^ VA 

T 3h \VA vp 



(55) 



(56) 



is the universal-limit thermal conductivity for a d-wave 
superconductor (without charge order) and we have de- 
fined A = Tq (our parameter of disorder) and B = ip'^ 
(our parameter of charge order). Inserting our expres- 
sions for the e's and A's from Eq. (fTU]) and integrating 
over p, we can obtain the zero-temperature thermal con- 
ductivity as a function of ip, Fq, and a = vp/vA- 



IV. ANALYTICAL RESULTS: CLEAN, 
ISOTROPIC LIMIT 

In the clean (A = Fq ^ 0), isotropic (a — vp/vA — 1) 
limit, the integrals in Eq. (j52p can be performed analyt- 
ically, providing us with a closed-form expression for the 
thermal conductivity tensor as a function of the charge 
density wave order parameter, ^. Selecting tpc (the value 
of ^ at which the nodes vanish) as our energy unit, and 
for a = 1, Eq. ^TU\i becomes 



ei = Pi + 1 Ai = p2 
£2 = J52 + 1 A2 = pi 

It is then useful to make a change of variables to 



91 

92 



Pi ~P2 

Pi + P2 + 1 



(57) 



(58) 



such that 



ei = [qi + 92 + l)/2 Ai = (92 - qi - l)/2 

£2 = (92 - 91 + l)/2 A2 = (qi +q2- l)/2 (59) 

Note that this change of variables has a Jacobian of 1/2, 
such that J cPp ^ 5 / (Pq- Therefore, 



(fqNi+ 7?,;iV2 

87r D 



(60) 
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iVi = AA 



A + B + 



where 



N2 = AAB [{q2 + If - ql] 



D= [f + A{q^ + 1 + 2B) + A^Y 



(61) 



(62) 



(63) 



(64) 



In the A — > limit, the numerator vanishes, so contribu- 
tions to the integral come only from the vicinity of points 
in q-space where the denominator vanishes as well, which 
requires / = 0. It is clear from Eq. that / is only 
equal to zero when q = 1 and q2 — B, the intersection 
of a unit circle about the origin and a horizonal line at 
q2 = B. 

For B > 1, there is no intersection, so the integral 
is zero. This is quite physical, since for i? > 1, ^ > 
ipc and the energy spectrum is gapped. Thus, in the 
clean, zero-temperature limit, there are no quasiparticles 
to transport heat and the thermal conductivity is zero. 

For _B < 1, the circle and line intersect at two points, 
q„ = (±-\/l — B^, B). These points are precisely the 
node and ghost-node of the energy spectrum, which will 
collide when ijj reaches tpc- For vanishing A, terms in A'^i, 
N2, and D that are higher than first order in A can be 
safely neglected and terms first order in A can be replaced 
by their values at q = q„. Doing so, we find that 



^ = {l + rj,B^)8{l + B)h 



where 



A 



8tt [f + 2A{l + B)Y 



(65) 



(66) 



and / is the function of q given in Eq. ((64|) . Changing 
variables to 



we see that 



xi — qi — 1 ~ X cos ( 
X2 = q2 — X sin 9 

?2 , 2 .3 . 



(67) 



/ = x'^/A + B^ + x-^ + x-^ COS0 -2Bxsme 
= ^ [1 + + 2/i(cos 9 cos 6*0 - sin 9 sin 6*0)] 

„2 



^ [l + h"^ + 2hcos{9 + eo)] 



(68) 



where 



h = 



2B 

and tan^o = (69) 



Then plugging / into Eq. ((66| . shifting 9 9 - 9o + n, 
and defining 7 = 2(1 -I- B)h^/x'^, wc find that 



1 f°° r h'^ 
Ii = — / dxx I d9 



A 



^ a;4 [i + /i2_ 2/1 cos 61 + ^7]' 



where 



87r(l + 5)7(5 ^ 



d9- 



A^l 



[1 + + yl7 - 2/icos( 



(70) 



(71) 



This integral over 9 is standard and easily evaluated via 
integration tablc^^. Doing so yields 



1 + ^^ 

^2=27r^^^i?(/i-l,^7) 



where 



D{u,T) 



rV2 



(w2+r2)3/2 



(72) 



(73) 



Since 7 is finite for all x, A7 vanishes as ^ ^ 0. There- 
fore, noting that 



lim D 

r->o 



(",r) = { 



for u^Q 
00 for u = 



and 



/oo 
duD{u,r) ^ 1 
-oc 



(74) 



(75) 



we see that D{u, F 0) is a representation of the Dirac 
delta function. Hence, 



and 



h = 



l2^-5{h-l) 



dx ■ 



16(1 + 5)70 x'^/A + B' 



(76) 



■5{h~l) (77) 



Noting that 5{h - 1) = 25{h'^ - 1) (since h > I) and 
letting u = a;2/2, this becomes 



1 



du ^ / 2u 



l + B)Jo u^ + B^ \u'^ + B'^ 



1 (78) 



For _B > 1, the argument of the delta function is never 
zero, so Ii — 0, as expected. For B < 1, the argument is 
zero at two points, u = u± = IzL ^/l — B^, so after a bit 
of delta- function gymnastics, we find that 



1 



1 



16(1 + B)y/T^rB^ Jo 

1 6(1-5) 
8(1 + 5) ^^-^W 



du [S{u — u_) + S{u — U-|_)] 
(79) 
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where Q{x) is the Heaviside step function. Finally, plug- 
ging back into Eq. ([55)) . we obtain a very simple, closed- 
form expression for the zero-temperature thermal con- 
ductivity tensor in the clean, isotropic limit. 



^0^ 
'^00 



1 + (V'/^c)^ 



e(^c-V') 







(80) 



(81) 



(82) 



These results are plotted in Fig. 21 For ip — 0, we re- 
cover the universal- limit thermal conductivity of a d-wave 
superconductor^ (see Eq. ([55]) V And for tp > tpc, as ex- 
pected, the thermal conductivity vanishes since the sys- 
tem has become gapped and there are no quasiparticles to 
transport the heat. For "0 between zero and ipc, thermal 
transport in the x and y directions differ, which makes 
sense as square symmetry has been explicitly broken by 
the charge density wave oriented in the x-direction. Par- 
allel to the CDW wave vector, thermal conductivity in 
the a;-direction decreases monotonically with ip, vanish- 
ing continuously at V'c- Perpendicular to the CDW wave 
vector, thermal conductivity in the y-direction increases 
with ip, exhibiting a square-root divergence before van- 
ishing abruptly at ipc- This divergence, a consequence 
of the clean limit, is replaced by a peak in Kq^ when 
nonzero disorder is considered, as will be shown in the 
next section. 



y° 2 



i 1 1 i 1 1 1 i i 




K / 
yy 




XX \ 








0.0 



0.5 1.0 



1.5 



FIG. 4: Calculated zero-temperature thermal conductivity 
tensor in the clean (Fo — > 0), isotropic {vf ~ va) limit. We 
plot and At™ as functions of the charge density wave order 
parameter, ip, from the closed-form expressions in Eqs. (I80p 
and (|81|) . As i/' approaches ^pc, the value beyond which the 
quasiparticle spectrum becomes gapped, 
uously while Kq^ diverges before dropping to zero. 



to, but distinct from, that of disorder, further smoothing 
the transition to zero thermal conductivity that occurs 
abruptly dA, ip = ipc in the clean, isotropic case. 



VI. CONCLUSIONS 



V. NUMERICAL RESULTS 

For the general case of nonzero disorder (Fo 7^ 0) 
and/or anisotropic Dirac nodes (a = vp/v\ 7^ 1), the 
p-space integration in Eq. ([5^ is more complicated, but 
can be computed numerically. Doing so, we calculated 
the zero-temperature thermal conductivity tensor as a 
function of charge density wave order parameter, ip, our 
parameter of disorder, Fg, and the anisotropy of the Dirac 
nodes, a = vp/vA- Results are plotted in Figs. [S] and [HI 

Fig. [S] shows Kq^ and as functions of ip for several 
values of Fq and a = 1. The clean-limit results calculated 
in Sec.|lV]are included (solid hues) for comparison. Note 
that as disorder increases, the transition to zero ther- 
mal conductivity at the nodal collision point (ip = ipc) 
gets rounded out, and the peak in Kq^ just prior to the 
collision point is diminished and broadened. Essentially, 
and not unexpectedly, disorder blurs the nodal collision, 
smoothing out the sharp transition seen in the clean case. 

The Fo = O.Obtpc results are reproduced in Fig. [SI along 
with plots of Kq^ and versus tp for larger values of a. 
For constant disorder, increasing a changes the shape of 
the Kq^ curve and diminishes and broadens the peak in 
Kq^ . The effect of increased nodal anisotropy is similar 



The coexistence of rf-wave superconductivity with 
charge order of sufficient magnitude can have a signif- 
icant effect on the energy spectrum of the Bogoliubov 
quasiparticles and the transport of heat by those quasi- 
particles at low temperatures. In this paper, we have 
considered a particularly simple form of charge order, a 
conventional s-wave charge density wave of wave vector 
Q = (tt/o, 0), the magnitude of which is characterized 
by a real, fc-independent order parameter, ip. The charge 
order halves the Brillouin zone, and as a function of 1/', 
the four nodes of the quasiparticle energy spectrum move 
in fc-space, approaching the reduced Brillouin zone edge. 
When ip reaches ipc, equal to the Fermi velocity times 
the fc-space distance from the original node location to 
the (7r/2,7r/2) point, the nodes reach the reduced Bril- 
louin zone edge and collide with their counterparts in 
the second reduced Brillouin zone. Beyond this point, 
the nodes vanish and the quasiparticle energy spectrum 
is fully gapped. 

We have used a linear response Kubo formula approach 
to calculate the zero temperature limit of the thermal 
conductivity tensor for this system. Working within an 
extended-Nambu basis (particle, hole, particle shifted by 
Q, hole shifted by Q), we constructed a 4 x 4 matrix 
Hamiltonian, Green's function, and thermal current op- 



9 



1.0 



0.8 



o 0.6 - 



0.4 



0.2 



0.0 - 



— r^=0 (Clean) 

• r„/v^=0.05 
□ r„/vi/=0.15 

* r /vi) =0.25 






0.5 1.0 




FIG. 5: Disorder dependence of calculated zero-temperature 
thermal conductivity tensor. We plot Kq^ (upper panel) and 



FIG. 6: Nodal anisotropy dependence of calculated zero- 
temperature thermal conductivity tensor. For fixed disorder 



(lower panel) as functions of charge density wave order (jpo = 0.05tpc), we plot Kg" (upper panel) and (lower 



parameter, tp, for several values of the disorder parameter, To 
In all cases, a — vf/va ~ 1. Included for comparison is the 
clean limit result (solid lines). Note that disorder smoothes 
the transition to zero thermal conductivity that results from 
the gapping of the energy spectrum a,t tp — tpc- The diver- 
gence of seen in the clean case is replaced by a peak that 
is diminished and broadened with increasing disorder. 



erator. We then used the Matsubara technique to eval- 
uate the bare-bubble thermal current-current correlator, 
neglecting vertex corrections and including disorder in 
the self-energy via a single broadening parameter, Tq. 
From this we calculated k^^ /T and k^^/T, in the limit 
of zero temperature, as a function of ip, Tq, and the nodal 
anisotropy a — vf/va- 

In the clean (Tq — s- 0), isotropic (vp = va) limit, our 
calculations yield a closed-form solution for the thermal 
conductivity tensor (plotted in Fig. 2]) 



(83) 



panel) as functions of charge density wave order parameter, 
'0, for several values of a = vf /va- Lines connecting the data 
points are guides to the eye. Note that the nodal transition at 
^ = t/'c is smoothed out by increasing velocity anisotropy in 
a manner similar to, but distinct from, the effect of disorder. 



^yy 



1 + (^/^c)^ 
Vl-(V'/V'c)^ 



where 



^00 _ k% IVF ^ VA 

T 'ih\vA vf 



(84) 



(85) 



(86) 



is the zero-temperature thermal conductivity for a d-wave 
superconductor with no charge order. As expected, the 
thermal conductivity takes the pure d-wave supercon- 
ductor value for ^ = and drops to zero for tjj > tpc. 
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where the quasiparticle energy spectrum has become fully 
gapped. For intermediate values ofip, Kq^ and k^^ differ, 
as square symmetry has been broken by the charge den- 
sity wave. For transport in the direction of the charge 
density wave vector, vanishes continuously as ip ap- 
proaches ipc- By contrast, for transport perpendicular 
to the charge density wave vector, Kq^ diverges before 
dropping abruptly to zero at ipc- This divergence is a 
consequence of the clean limit and is replaced by a finite 
peak in the presence of disorder. 

For the more complicated case of nonzero disorder 
(To ^ 0) and/or anisotropic nodes {vp 7^ va), we have 
obtained results via a numerical calculation. We find 
that disorder smoothes out the transition to zero thermal 
conductivity across the nodal collision (see Fig. [5]). The 
clean-limit divergence in Kq^ just before the transition is 
replaced by a peak which broadens and decreases in am- 
plitude with increasing disorder. The abrupt drop in the 
clean-limit Kq^ is similarly broadened. Essentially, the 
disorder- broadening of the quasiparticle spectral function 
averages over what was, in the clean limit, a sharp transi- 
tion from gapless to gapped quasiparticles. We find that 
increased nodal anisotropy has a similar effect, amplify- 
ing the disorder effect and thereby further broadening the 
features in the ■(/'-dependence of the thermal conductiv- 
ity (see Fig. [5]). And the fact that disorder has an effect 
indicates that the low-temperature thermal conductivity 
is no longer universal (disorder-independent) in the pres- 
ence of charge order, which is in line with the results of 
recent measurement o^^'^^i'^^'^^ of low-temperature ther- 
mal transport in the underdoped cuprates. 

In these calculations, we have enjoyed the theorist's 
luxury of being able to turn on, by hand, a charge den- 
sity wave to coexist with the d-wave superconductiv- 
ity. The experimenter does not have direct access to 
such a knob. However, in the d-wave superconduct- 
ing state of the cuprates, charge order does appear to 
be enhanced with underdoping. Hence the features of 
the V'-dependent thermal conductivity curves calculated 
herein should serve as signatures for the underdoping- 
dependence of thermal conductivity measured in the un- 
derdoped cuprates. Of course, most dramatic would 
be the observation of the nodal collision beyond which 
the low-temperature thermal conductivity drops to zero. 
However, even if the amplitude of charge order is insuf- 
ficient to reach the nodal collision, these results should 
provide insight to the approach to the transition. 

A sequel to this work, exploring the effects of a more 
elaborate model of disorder, as well as the contribution of 
vertex corrections, is in preparation''"'^. Future work will 
also examine the effect of different types of charge order 
(beyond the conventional s-wave case considered here) 
of different wave vector (beyond the unit-cell-doubling 
Q = {n/a,0) case considered here) and of multiple wave 
vectors (like the checkerboard charge order observed in 
some cupratc G^^i^^ ). 
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APPENDIX: SUBTLETIES OF THE SPECTRAL 
REPRESENTATION 

In the calculations described in this paper, we have 
made use of the 4x4 extended-Nambu basis of Eq. ([5]) 
for the Hamiltonian and Green's functions. This choice of 
basis provides a compact realization of the Hamiltonian 
and is quite convenient in many respects. However, use of 
a matrix Green's function does introduce some subtleties 
regarding the spectral representation, and we would like 
to address those here. 

All of our results could have been obtained by diag- 
onalizing the Hamiltonian from the outset and working 
with the diagonalized Green's function 



Gz5(itj) = U''G(iuj)U 
with diagonal matrix elements 



[Gd(*c.)]„„ = 



1 



ILU — E 



(n) 



E(«)(iw) 



(A.l) 



(A.2) 



where the i?^." are the eigenvalues of Hk, the E*^") are 
the corresponding self-energies, and the eigenvectors de- 
fine the columns of unitary transformation matrix U. In 
the diagonal basis, it is quite valid to define a spectral 
representation for the Green's function 



duj 



' iuj — ijJi 



(A.3) 



where 



AdH EE ^ [Gg(c.) - G^u;)] = -ilmGg(c.) (A.4) 

Note that the second equality follows from the fact 
that the retarded diagonal Green's function, G^{uj) = 
Goii^ ~* oj -\- is the complex conjugate of the ad- 
vanced diagonal Green's function, G^{lo) = Gnii^ 
Lo — i5), which is clear from Eq. (|A.2[) . 

The non-diagonal matrix Green's function can there- 
fore be expressed as 



G{iuj) = UGD{iuj)U^ = / duji 



-iC/ImGg(c^i)C/t 



IbJ — LUl 



(A.5) 



which is not equivalent to the right-hand side of Eq. 



^^^-llmG^(.0 r^^^-llm[^Gg(.O^t] 



tUJ — U!i 



lOJ — UJi 



(A.6) 
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unless the diagonalization transformation commutes with 
taking the imaginary part. Equivalently, note that it is 
vahd to define a non-diagonal matrix spectral function, 

A{uj) = UAd{uj)W, such that 



G{iuj) = 




(A.7) 



but 

AH^_L[G«(^)„G^(^)] (A.8) 



will not be equal to — ImG^/Tr (or even have to be 
real) unless the non-diagonal retarded Green's function, 
G^{uj) = G{iLjj Lo + iS), is the complex conjugate of 
the non-diagonal advanced Green's function, G"^{uj) = 
G(iuj ^ oj — iS). 

For the case of real ip that we consider in this paper, 
Hk and therefore U are real, so diagonalization does com- 
mute with taking the imaginary part, G^ is the complex 
conjugate of G^, and Eq. is valid. But this is not 
generically the case for complex tp. 
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